Statistical analysis

baRulho:quantifying habitat-induced degradation of (animal) acoustic signals
University of Costa Rica

Author

Marcelo Araya-Salas

Published

¶ “2023-04-20” ¶

Source code, data and annotation protocol found at

1 Purposes

  • Estimate degradation metrics on re-recorded signals from playback experiment at Bosque de Tlalpan, Mexico City, 2019

  • Run statistical analyses

 

Load packages

Code 
# install/ load packages
sketchy::load_packages(packages = c("remotes", "kableExtra", "knitr",
    "formatR", "rprojroot", "maRce10/warbleR", "ggplot2", "tidyr",
    "viridis", "corrplot", "brms", "ggdist", "cowplot", "cmdstanr",
    "maRce10/brmsish"))

source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
Code 
degrad_df <- read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")

degrad_params <- c("blur.ratio", "spectrum.blur.ratio", "envelope.correlation",
    "excess.attenuation", "signal.to.noise.ratio", "cross.correlation",
    "tail.to.signal.ratio", "tail.to.noise.ratio", "spectrum.correlation")

comp.cases <- complete.cases(degrad_df[, names(degrad_df) %in% degrad_params])

pca <- prcomp(degrad_df[comp.cases, names(degrad_df) %in% degrad_params],
    scale. = TRUE)

# add to data
degrad_df$pc1 <- NA
degrad_df$pc1[comp.cases] <- pca$x[, 1]
degrad_df$pc1.1m.rate <- degrad_df$distance
# plot rotation values by PC
pca_rot <- as.data.frame(pca$rotation[, 1:4])

pca_rot_stck <- stack(pca_rot)

pca_rot_stck$variable <- rownames(pca_rot)
pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative")
pca_rot_stck$rotation <- abs(pca_rot_stck$values)

ggplot(pca_rot_stck, aes(x = variable, y = rotation, fill = Sign)) +
    geom_col() + coord_flip() + scale_fill_viridis_d(alpha = 0.7,
    begin = 0.3, end = 0.8) + facet_wrap(~ind)

1.1 Correlation among metrics

Raw metrics:

Code 
cormat <- cor(degrad_df[, degrad_params], use = "pairwise.complete.obs")

rownames(cormat) <- colnames(cormat) <- degrad_params

cols_corr <- colorRampPalette(c("white", "white", viridis(4, direction = -1)))(10)

cp <- corrplot.mixed(cormat, tl.cex = 0.7, upper.col = cols_corr,
    lower.col = cols_corr, order = "hclust", lower = "number", upper = "ellipse",
    tl.col = "black")

Code 
# sort parameters as in clusters for cross correlation
degrad_params <- degrad_params[match(rownames(cp$corr), degrad_params)]

1.2 Data description

  • 20 frequencies
  • 3 locations
  • 11520 test sounds
  • 160 sound treatment combinations
  • Sample sizes per location, transect and signal type
Code 
agg <- aggregate(cbind(sound.id, treatment.replicates) ~ location +
    habitat.structure + distance, degrad_df, function(x) length(unique(x)))

agg$replicates <- agg$sound.id/agg$treatment.replicates

agg
location habitat.structure distance sound.id treatment.replicates replicates
1 closed 10 480 160 3
2 closed 10 480 160 3
3 closed 10 480 160 3
1 open 10 480 160 3
2 open 10 480 160 3
3 open 10 480 160 3
1 closed 30 480 160 3
2 closed 30 480 160 3
3 closed 30 480 160 3
1 open 30 480 160 3
2 open 30 480 160 3
3 open 30 480 160 3
1 closed 65 480 160 3
2 closed 65 480 160 3
3 closed 65 480 160 3
1 open 65 480 160 3
2 open 65 480 160 3
3 open 65 480 160 3
1 closed 100 480 160 3
2 closed 100 480 160 3
3 closed 100 480 160 3
1 open 100 480 160 3
2 open 100 480 160 3
3 open 100 480 160 3

1.3 Statistical analysis

Code 
iter <- 10000
chains <- 4
priors <- c(prior(normal(0, 6), class = "b"), prior(cauchy(0, 4),
    class = "sd"))


# set frequency to mean-centered
degrad_df$frequency <- scale(degrad_df$frequency, scale = FALSE)

# set base level for factors
degrad_df$habitat.structure <- factor(degrad_df$habitat.structure,
    levels = c("open", "closed"))
degrad_df$frequency.modulation <- factor(degrad_df$frequency.modulation,
    levels = c("no_fm", "fm"))
degrad_df$amplitude.modulation <- factor(degrad_df$amplitude.modulation,
    levels = c("no_am", "am"))
degrad_df$duration <- factor(degrad_df$duration, levels = c("short",
    "long"))
degrad_df$location <- as.factor(degrad_df$location)
# degrad_df$ord.distance <- ordered(degrad_df$distance)

set.seed(123)

cmdstanr::set_cmdstan_path("~/Documentos/cmdstan/")

# to run within-chain parallelization
mod_pc1 <- brm(formula = pc1 ~ frequency * habitat.structure + frequency.modulation *
    habitat.structure + amplitude.modulation * habitat.structure +
    duration * habitat.structure + mo(distance) + (1 | location) +
    (1 | treatment.replicates), data = degrad_df, prior = priors,
    iter = iter, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/regression_model_pc1.RDS",
    file_refit = "always", backend = "cmdstanr", threads = threading(10))

mod_blurratio <- brm(formula = blur.ratio ~ frequency * habitat.structure +
    frequency.modulation * habitat.structure + amplitude.modulation *
    habitat.structure + duration * habitat.structure + mo(distance) +
    (1 | location) + (1 | treatment.replicates), data = degrad_df,
    prior = priors, iter = iter, chains = chains, cores = chains,
    control = list(adapt_delta = 0.99, max_treedepth = 15), file = "./data/processed/regression_model_blur_ratio.RDS",
    file_refit = "always", backend = "cmdstanr", threads = threading(10))
Code 
my.viridis <- function(...) viridis(alpha = 0.5, begin = 0.3, end = 0.7,
    ...)


effects <- c("habitat structure", "frequency modulation", "amplitude modulation",
    "frequency", "duration", "frequency:habitat structure", "habitat structure:duration",
    "habitat structure:amplitude modulation", "habitat structure:frequency modulation")
extended_summary(read.file = "./data/processed/regression_model_pc1.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    )

1.4 regression_model_pc1

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, -0.3, 2.6) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.6) simo-dirichlet(1) pc1 ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance) + (1 | location) + (1 | treatment.replicates) 10000 4 1 5000 36 0 2413.73 5075.71 721775339
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
frequency 0.155 0.126 0.184 1.001 2714.855 5483.055
habitat structure 1.065 0.981 1.148 1 15427.177 14468.772
frequency modulation 0.809 0.643 0.975 1.002 2499.781 5266.845
amplitude modulation 0.233 0.066 0.401 1.002 2483.783 5075.710
duration 0.151 -0.018 0.318 1 2413.730 5470.687
frequency:habitat structure 0.055 0.041 0.070 1 28532.118 14111.309
habitat structure:frequency modulation -0.154 -0.238 -0.071 1 23138.819 14435.455
habitat structure:amplitude modulation -0.017 -0.100 0.066 1 23604.229 14414.408
habitat structure:duration 0.080 -0.003 0.164 1 21520.640 14884.870

Code 
extended_summary(read.file = "./data/processed/regression_model_blur_ratio.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects)

1.5 regression_model_blur_ratio

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, 0.1, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) blur.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance) + (1 | location) + (1 | treatment.replicates) 10000 4 1 5000 83 0 2401.559 4924.56 721775339
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 0.029 0.026 0.033 1 19677.629 14582.926
frequency modulation 0.065 0.059 0.072 1.002 2440.788 5134.879
amplitude modulation 0.024 0.017 0.031 1.002 2571.274 5338.939
frequency 0.001 0.000 0.002 1.001 3034.127 6220.818
duration 0.001 -0.006 0.008 1.002 2401.559 4924.560
frequency:habitat structure 0.003 0.002 0.003 1 20247.850 12329.382
habitat structure:duration 0.005 0.001 0.009 1 29331.488 15213.600
habitat structure:amplitude modulation 0.002 -0.001 0.006 1 30029.842 15156.579
habitat structure:frequency modulation -0.020 -0.024 -0.017 1 29311.245 15487.350

Code 
extended_summary(read.file = "./data/processed/regression_model_excess_attenuation.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects)

1.6 regression_model_excess_attenuation

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, 83.8, 34) sd-cauchy(0, 4) sigma-student_t(3, 0, 34) simo-dirichlet(1) excess.attenuation ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance) + (1 | location) + (1 | treatment.replicates) 10000 4 1 5000 3 0 1918.559 3646.097 497531395
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 9.887 9.056 10.729 1 11112.681 13364.153
frequency modulation -3.197 -4.695 -1.673 1.002 1918.559 3646.097
amplitude modulation -2.676 -4.217 -1.174 1.001 2235.177 3941.118
frequency 3.163 2.902 3.426 1.001 2051.986 4589.126
duration -0.896 -2.424 0.606 1.001 2053.495 4359.213
frequency:habitat structure -0.639 -0.786 -0.493 1 20178.460 15154.954
habitat structure:duration -0.391 -1.231 0.446 1 16298.977 14986.898
habitat structure:amplitude modulation -1.575 -2.406 -0.746 1 16186.223 14732.686
habitat structure:frequency modulation -0.127 -0.957 0.711 1 16971.731 14743.070

Code 
extended_summary(read.file = "./data/processed/regression_model_signal-to-noise_ratio.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects)

1.7 regression_model_signal-to-noise_ratio

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, 5.4, 7.6) sd-cauchy(0, 4) sigma-student_t(3, 0, 7.6) simo-dirichlet(1) signal.to.noise.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance) + (1 | location) + (1 | treatment.replicates) 10000 4 1 5000 11 0 1662.527 3288.461 865228054
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure -2.048 -2.342 -1.755 1 8927.894 11843.910
frequency modulation 0.975 0.374 1.573 1.002 1726.048 3288.461
amplitude modulation -1.229 -1.835 -0.620 1.002 1745.506 3480.490
frequency -0.441 -0.546 -0.339 1.001 1662.527 3899.319
duration 0.917 0.322 1.514 1.002 1715.016 3422.535
frequency:habitat structure 0.099 0.048 0.149 1 19636.322 14899.590
habitat structure:duration -0.189 -0.476 0.107 1 14528.607 14668.777
habitat structure:amplitude modulation 0.481 0.191 0.773 1 14124.868 14296.795
habitat structure:frequency modulation -0.669 -0.964 -0.372 1 13248.444 14006.836

Code 
extended_summary(read.file = "./data/processed/regression_model_tail-to-signal_ratio.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects)

1.8 regression_model_tail-to-signal_ratio

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, -6.3, 8.7) sd-cauchy(0, 4) sigma-student_t(3, 0, 8.7) simo-dirichlet(1) tail.to.signal.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance) + (1 | location) + (1 | treatment.replicates) 10000 4 1 5000 3 0 2040.207 4477.209 1732769608
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 3.253 2.898 3.608 1 10995.132 12205.423
frequency modulation -0.817 -1.434 -0.207 1.001 2040.207 4477.209
amplitude modulation 0.814 0.207 1.435 1.002 2255.951 4695.752
frequency 0.302 0.196 0.413 1.002 2322.046 4559.236
duration -0.427 -1.049 0.183 1.003 2433.154 5441.122
frequency:habitat structure 0.020 -0.041 0.082 1.001 20915.048 14502.400
habitat structure:duration -0.055 -0.408 0.296 1 16428.408 14629.440
habitat structure:amplitude modulation -0.386 -0.740 -0.032 1 16037.645 14744.846
habitat structure:frequency modulation 0.542 0.187 0.901 1 15537.887 14213.645

 

2 Takeaways

 

 

Session information

R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] brmsish_1.0.0      cmdstanr_0.4.0     cowplot_1.1.1      ggdist_3.2.0      
 [5] brms_2.18.0        Rcpp_1.0.9         corrplot_0.90      viridis_0.6.2     
 [9] viridisLite_0.4.1  tidyr_1.1.3        ggplot2_3.4.0      warbleR_1.1.28    
[13] NatureSounds_1.0.4 seewave_2.2.0      tuneR_1.4.1        rprojroot_2.0.3   
[17] formatR_1.11       knitr_1.42         kableExtra_1.3.4   remotes_2.4.2     

loaded via a namespace (and not attached):
  [1] uuid_1.1-0           backports_1.4.1      systemfonts_1.0.4   
  [4] plyr_1.8.7           igraph_1.3.5         splines_4.1.0       
  [7] crosstalk_1.2.0      TH.data_1.1-0        rstantools_2.2.0    
 [10] inline_0.3.19        digest_0.6.31        htmltools_0.5.4     
 [13] fansi_1.0.3          magrittr_2.0.3       checkmate_2.1.0     
 [16] RcppParallel_5.1.5   matrixStats_0.62.0   xts_0.12.2          
 [19] sandwich_3.0-1       svglite_2.1.0        prettyunits_1.1.1   
 [22] colorspace_2.0-3     signal_0.7-7         rvest_1.0.3         
 [25] xfun_0.36            dplyr_1.0.10         callr_3.7.3         
 [28] crayon_1.5.2         RCurl_1.98-1.9       jsonlite_1.8.4      
 [31] lme4_1.1-27.1        ape_5.6-2            survival_3.2-11     
 [34] zoo_1.8-11           glue_1.6.2           gtable_0.3.1        
 [37] emmeans_1.8.1-1      webshot_0.5.4        distributional_0.3.1
 [40] pkgbuild_1.4.0       rstan_2.21.7         abind_1.4-5         
 [43] scales_1.2.1         mvtnorm_1.1-3        DBI_1.1.1           
 [46] xaringanExtra_0.7.0  miniUI_0.1.1.1       dtw_1.23-1          
 [49] xtable_1.8-4         proxy_0.4-27         stats4_4.1.0        
 [52] StanHeaders_2.21.0-7 DT_0.26              htmlwidgets_1.5.4   
 [55] httr_1.4.4           threejs_0.3.3        posterior_1.3.1     
 [58] ellipsis_0.3.2       pkgconfig_2.0.3      loo_2.4.1.9000      
 [61] farver_2.1.1         utf8_1.2.2           labeling_0.4.2      
 [64] tidyselect_1.2.0     rlang_1.0.6          reshape2_1.4.4      
 [67] later_1.3.0          munsell_0.5.0        tools_4.1.0         
 [70] cli_3.6.0            generics_0.1.3       ggridges_0.5.4      
 [73] evaluate_0.20        stringr_1.5.0        fastmap_1.1.0       
 [76] yaml_2.3.7           processx_3.8.0       purrr_1.0.0         
 [79] packrat_0.9.0        pbapply_1.6-0        nlme_3.1-152        
 [82] projpred_2.0.2       mime_0.12            xml2_1.3.3          
 [85] compiler_4.1.0       bayesplot_1.9.0      shinythemes_1.2.0   
 [88] rstudioapi_0.14      gamm4_0.2-6          sketchy_1.0.2       
 [91] tibble_3.1.8         stringi_1.7.12       highr_0.10          
 [94] ps_1.7.2             Brobdingnag_1.2-9    lattice_0.20-44     
 [97] Matrix_1.5-1         nloptr_1.2.2.2       markdown_1.3        
[100] shinyjs_2.1.0        fftw_1.0-7           tensorA_0.36.2      
[103] vctrs_0.5.2          pillar_1.8.1         lifecycle_1.0.3     
[106] bridgesampling_1.1-2 estimability_1.4.1   bitops_1.0-7        
[109] httpuv_1.6.6         R6_2.5.1             promises_1.2.0.1    
[112] gridExtra_2.3        codetools_0.2-18     boot_1.3-28         
[115] colourpicker_1.2.0   MASS_7.3-54          gtools_3.9.3        
[118] assertthat_0.2.1     rjson_0.2.21         withr_2.5.0         
[121] shinystan_2.6.0      multcomp_1.4-17      mgcv_1.8-36         
[124] parallel_4.1.0       grid_4.1.0           minqa_1.2.4         
[127] coda_0.19-4          rmarkdown_2.20       shiny_1.7.3         
[130] base64enc_0.1-3      dygraphs_1.1.1.6